Methods and systems for quantifying coherency and constraining coherency-based separation in simultaneous shooting acquisition

ABSTRACT

This disclosure presents methods and systems for deblending blended seismic data obtained during simultaneous shooting acquisition into deblended seismic data gathers. Methods and systems iteratively separate the blended seismic data into the deblended seismic data gathers based on semblance analysis of a residual difference between the blended seismic data and the deblended seismic data gathers. Each deblended seismic data gather is associated with one of the sources and appears to have been obtained without substantial interference from seismic energy produced by other sources.

CROSS-REFERENCE TO RELATED APPLICATION

This application claims the benefit of provisional application No.61/944,688, filed Feb. 26, 2014.

BACKGROUND

In the past few decades, the petroleum industry has invested heavily inthe development of marine seismic survey techniques that yield knowledgeof subterranean formations beneath a body of water in order to find andextract valuable petroleum resources. A typical marine seismic surveymay be carried-out with a survey vessel towing a seismic source and thesame vessel, or another vessel, towing one or more streamers that form aseismic data acquisition array below the surface of the water and abovethe subterranean formation. The survey vessel typically contains seismicdata acquisition equipment, such as navigation control, seismic sourcecontrol, seismic receiver control, and data recording equipment. Theseismic source control activates a seismic source, which is typically anarray of source elements, such as air guns or marine vibrators, thatproduces acoustic signals at selected times or at selected locations.Each acoustic signal is generally a sound wave that travels down throughthe body of water and into the subterranean formation. At interfacesbetween different types of rock, a portion of the sound wave istransmitted, a portion is refracted, and a portion is reflected backinto the body of water to propagate toward the water surface. Thestreamers towed behind the survey vessel are elongated cable-likestructures. Each streamer may include a number of seismic receivers ormulti-component sensors that detect pressure and/or particle motionwavefields of the sound waves reflected back into the water from thesubterranean formation. However, the wavefields measured by thereceivers may also include seismic energy created by other activesources. This additional seismic energy may create seismic interference.Those working in the petroleum industry seek techniques to improveseparation of seismic energy created by two or more sources.

DESCRIPTION OF THE DRAWINGS

FIGS. 1A-1B show side-elevation and top views, respectively, of anexample data acquisition system.

FIG. 2 shows an example of simultaneous shooting data acquisition.

FIG. 3 shows a side-elevation view of a data acquisition system with amagnified view of a receiver.

FIG. 4 shows example rays that represent acoustic signal paths outputfrom a source.

FIG. 5 shows an example shot gather with five traces.

FIG. 6 shows an example of a gather with thirty-eight traces.

FIG. 7 shows a plot of different ways seismic data may be sorted intodomains.

FIG. 8 shows a mathematical representation of a relationship betweenblended and deblended seismic data.

FIG. 9 shows an example mathematical representation of a relationshipbetween blended and deblended seismic data obtained using the exampledata acquisition system shown in FIG. 1.

FIG. 10 shows an example of a deblended data matrix.

FIG. 11 shows an example flow-control diagram of an iterative method forseparating seismic data obtained using simultaneous shootingacquisition.

FIG. 12 shows a flow-control diagram for a routine “update residual anddeblended seismic data” called in the flow-control diagram of FIG. 11.

FIG. 13 shows an example of a generalized computer system.

FIGS. 14A-14C show results of seismic source separation applied forcommon-shot gathers.

DETAILED DESCRIPTION

In simultaneous shooting acquisition (“SSA”), two or more sources areactivated in succession with dithered time delays above a subterraneanformation. The term “simultaneous” does not mean that the two or moresources are activated at exactly the same time. SSA refers to activationof the two or more sources with different time delays in a short timeinterval. As used herein, “simultaneous” shooting or activation meansactivating two or more sources within a short time interval that, forexample, is small compared to the seismic recording period. Activationof the two or more sources within a short time interval may createinterfering seismic energy emitted from the subterranean formation. Adata acquisition system may record the interfering seismic energy asblended seismic data. One aim of SSA is to reduce the time to acquireseismic data and/or to increase the diversity of the seismic data interms of fold, azimuth, and offsets. This disclosure presents methodsand systems that may be used to deblend the blended seismic data intoseparate deblended seismic data gathers. Each deblended seismic datagather may be associated with one of the activated sources, and mayappear to have been obtained without substantial interference fromseismic energy produced by other sources.

Methods of the current disclosure iteratively separate the blendedseismic data into deblended seismic data gathers over the individualsources. In some implementations, significant coherent energy content inthe blended seismic data may be identified and separated. Once thesignificant coherent energy components have been identified andseparated, coherent energy events that are strongly interfered with byincoherent energy may become better accessible and may be separated moreeasily in subsequent iterations. A coherency filter, such as a medianfilter, frequency-space deconvolution, tau-p filter, and afrequency-wavenumber filter, may be used to filter coherent energy ineach iteration. Coherency filters are typically efficient when thecoherent energy is greater than the interfering incoherent energy.Typical separation techniques may leave artifacts, such as energyleakage or smearing effects, that contaminate deblended seismic datagathers when these same coherency filters are applied to parts ofblended seismic data with interfering incoherent energy that is greaterthan the coherent energy. On the other hand, methods of the currentdisclosure may produce deblended seismic data gathers that aresubstantially free of energy leakage or smearing effects by selectivelyapplying the coherency filters to parts of the seismic data wherecoherent energy is greater than incoherent energy, and leave unfilteredother parts of the same seismic data identified as having interferingincoherent energy that is greater than the coherent energy. Insubsequent iterations, when the interfering incoherent energy has beenseparated, the weaker coherent energy may be filtered/separated usingthe coherency filters.

Methods of the current disclosure use semblance analysis to bettercontrol the iterative separation process. Semblance analysis may be usedto generate a semblance value for each sample of seismic data. Thesemblance value may give a level of coherent energy for the sample ofseismic data and may be evaluated based on a threshold to determinewhether or not a coherency filter should be applied to the sample ofseismic data. For the first iteration, a semblance value of a sample ofseismic data may be small because of strong interfering seismic energy,while in subsequent iterations, when much of the interfering energy hasbeen reduced, the semblance value for the same sample may be increased.The threshold may be set as a user-defined process parameter. In certainimplementations, one threshold may be used for the full iterativeseparation process in which the separation process stops when no moreseismic data with a larger semblance value than the threshold isavailable. In other implementations, the threshold may be automaticallyadjusted as the number of iterations increases. In particular, a highthreshold may be set at the beginning of the separation process in orderto initially focus on the greatest coherent energy events in the seismicdata followed by progressively decreasing the threshold as the number ofiterations increases until the recorded seismic data is separated intothe separate deblended seismic data gathers.

FIGS. 1A-1B show side-elevation and top views, respectively, of anexample data acquisition system 100 that may be used for SSA. In thisexample, the data acquisition system is composed of four survey vessels101-104 towing six sources denoted by S1, S2, S3, S4, S5, and S6. Surveyvessel 101 also tows six separate streamers 106-111 beneath a freesurface 112 of a body of water. The body of water can be, for example,an ocean, a sea, a lake, or a river, or any portion thereof. In thisexample, each streamer is attached at one end to the survey vessel 101via a streamer-data-transmission cable. The illustrated streamers106-111 form a planar horizontal data acquisition surface with respectto the free surface 112. It should be noted that the number of sourcesused in SSA is not limited to six sources. In practice, the number ofsources S used in SSA may range from as few as two sources to more thansix sources, such as ten or more sources, located at various positionsin relation to the data acquisition surface. The data acquisitionsurface may also be smoothly varying due to active sea currents andweather conditions. In other words, although the streamers 106-111 areillustrated in FIGS. 1A and 1B and subsequent figures as straight andsubstantially parallel to the free surface 112, in practice, the towedstreamers may undulate as a result of dynamic conditions of the body ofwater in which the streamers are submerged. A data acquisition surfaceis not limited to having a planar horizontal orientation with respect tothe free surface 112. The streamers may be towed at depths that anglethe data acquisition surface with respect to the free surface 112 or oneor more of the streamers may be towed at different depths. A dataacquisition surface is not limited to six streamers as shown in FIG. 1B.In practice, the number of streamers used to form a data acquisitionsurface can range from as few as one streamer to as many as 20 or morestreamers. A marine seismic survey may also be conducted with additionalsurvey vessels (not show) towing additional data acquisition surfaceseach with up to 20 or more streamers. For example, the additional dataacquisition surfaces may be towed substantially parallel to the dataacquisition surface towed by the survey vessel 101.

FIG. 1A includes an xz-plane 114, and FIG. 1B includes an xy-plane 116of the same Cartesian coordinate system having three orthogonal, spatialcoordinate axes labeled x, y and z. The coordinate system is used tospecify orientations and coordinate locations within the body of water.The x-direction specifies the position of a point in a directionparallel to the length of the streamers (or a specified portion thereofwhen the length of the streamers are curved) and is referred to as the“in-line” direction. The y-direction specifies the position of a pointin a direction perpendicular to the x-axis and substantially parallel tothe free surface 112 and is referred to as the “cross-line” direction.The z-direction specifies the position of a point perpendicular to thexy-plane (i.e., perpendicular to the free surface 112) with the positivez-direction pointing downward away from the free surface 112. Thestreamers 106-111 are generally long cables containing power anddata-transmission lines that connect receivers 118 (represented byshaded rectangles) spaced-apart along the length of each streamer toseismic data acquisition equipment and data-storage devices that may belocated on board the survey vessel 101.

Streamer depth below the free surface 112 may be estimated at variouslocations along the streamers using depth-measuring devices attached tothe streamers. For example, the depth-measuring devices can measurehydrostatic pressure or utilize acoustic distance measurements. Thedepth-measuring devices can be integrated with depth controllers, suchas paravanes or water kites that control and maintain the depth andposition of the streamers as the streamers are towed through the body ofwater. The depth-measuring devices are typically placed at intervals(e.g., about 300 meter intervals in some implementations) along eachstreamer. Note that in other implementations buoys may be attached tothe streamers and used to maintain the orientation and depth of thestreamers below the free surface 112.

FIG. 1A shows a cross-sectional view of the survey vessels 101-104towing the sources S1-S6 above a subterranean formation 120. Curve 122,the formation surface, represents a top surface of the subterraneanformation 120 located at the bottom of the body of water. Thesubterranean formation 120 may be composed of a number of subterraneanlayers of sediment and rock. Curves 124, 126, and 128 representinterfaces between subterranean layers of different compositions. Ashaded region 130, bounded at the top by a curve 132 and at the bottomby a curve 134, represents a subterranean hydrocarbon deposit, the depthand positional coordinates of which may be determined, at least in part,by analysis of seismic data collected during a marine seismic survey.Each of the sources S1-S6 may be an air gun, marine vibrator, orcomposed of an array of air guns and/or marine vibrators. FIG. 1Aillustrates an acoustic signal expanding outward from the source S1 as apressure wavefield 136 represented by semicircles of increasing radiuscentered at the source S1. The outwardly expanding wavefronts from eachof the sources may be spherical but are shown in vertical plane crosssection in FIG. 1A. The outward and downward expanding portion of thepressure wavefield 136 is called the “primary wavefield,” whicheventually reaches the formation surface 122 of the subterraneanformation 120, at which point the primary wavefield may be partiallyreflected from the formation surface 122 and partially refracteddownward into the subterranean formation 120, becoming elastic waveswithin the subterranean formation 120. In other words, in the body ofwater, the acoustic signal is composed primarily of compressionalpressure waves, or P-waves, while in the subterranean formation 120, thewaves include both P-waves and transverse waves, or S-waves. Within thesubterranean formation 120, at each interface between different types ofmaterials or at discontinuities in density or in one or more of variousother physical characteristics or parameters, downward propagating wavesmay be partially reflected and partially refracted. As a result, eachpoint of the formation surface 122 and each point of the interfaces 124,126, and 128 may be a reflector that becomes a potential secondary pointsource from which acoustic and elastic wave energy, respectively, mayemanate upward toward the receivers 118 in response to the acousticsignal generated by the source S1 and downward-propagating elastic wavesgenerated from the pressure impulse. As shown in FIG. 1A, secondarywaves of significant amplitude may be generally emitted from points onor close to the formation surface 122, such as point 138, and frompoints on or very close to interfaces in the subterranean formation 120,such as points 140 and 142.

The secondary waves may be generally emitted at different times within arange of times following the initial acoustic signal. A point on theformation surface 122, such as the point 138, may receive a pressuredisturbance from the primary wavefield more quickly than a point withinthe subterranean formation 120, such as points 140 and 142. Similarly, apoint on the formation surface 122 directly beneath the source S1 mayreceive the pressure disturbance sooner than a more distant-lying pointon the formation surface 122. Thus, the times at which secondary andhigher-order waves are emitted from various points within thesubterranean formation 120 may be related to the distance, inthree-dimensional space, of the points from the activated source.

Acoustic and elastic waves, however, may travel at different velocitieswithin different materials as well as within the same material underdifferent pressures. Therefore, the travel times of the primarywavefield and secondary wavefield emitted in response to the primarywavefield may be functions of distance from each source as well as thematerials and physical characteristics of the materials through whichthe wavefields travel. In addition, the secondary expanding wavefrontsmay be altered as the wavefronts cross interfaces and as the velocity ofsound varies in the media are traversed by the wave. The superpositionof waves emitted from within the subterranean formation 120 in responseto the primary wavefield may be a generally complicated wavefield thatincludes information about the shapes, sizes, and materialcharacteristics of the subterranean formation 120, including informationabout the shapes, sizes, and locations of the various reflectingfeatures within the subterranean formation 120 of interest toexploration seismologists.

FIG. 2 shows an example of carrying out SSA using the data acquisitionsystem 100 traveling a ship track represented by a directional arrow202. Two or more of the sources S1-S6 may be activated at locationsrepresented by filled circles 204-206 as the data acquisition system 100travels along the ship track 202. Directional arrow 208 represents amarine-seismic-survey-timeline that indicates when the sources S1-S6 areactivated and when seismic data is recorded as the data acquisitionsystem 100 travels the ship track 202. The example timeline 208 issegmented into recording periods 210-212. At the beginning of each ofthe recording periods 210-212, two or more of the sources S1-S6 areactivated in shot intervals represented by boxes 214-216. For example,the shot intervals may have durations of approximately 1 second orlonger and the recording periods may range from approximately 6-20seconds. Seismic data may be recorded separately in each recordingperiod or recorded continuously as the data acquisition system 100travels the ship track 202. In other words, there may be no separationon the timeline 208 between two adjacent recording periods, or twoadjacent recording periods may be separated by a pause period (notshown). Pause periods may have durations from less than a second toabout 5 to 10 seconds.

In SSA, two or more of the sources S1-S6 are selected for activationwithin each shot interval. SSA refers to activation of the sources S1-S6at different times within the same short shot interval, which may resultin temporal overlap of seismic energy generated by each of the sources.As a result, the seismic data recorded in a recording period may be arecord of interfering seismic energy produced by each activated sourceand is called “blended seismic data.” Activating two or more of thesources S1-S6 in a single shot interval at randomly, pseudo-random, orpatterned (e.g., sinusoidal time delay) activation times is called“dithering.” Generally, two or more of the sources S1-S6 may beactivated with a random, pseudo-random, or patterned time delay, τ_(s),(i.e., dithered time delay) after the beginning of the shot interval,where subscript s is a positive integer source index. For example, FIG.2 includes example detailed views 218 and 220 of dithered sourceactivation times in the shot intervals 214 and 215, respectively. Theshot interval 214 begins when the source S1 is activated, whichcorresponds to a time delay of τ₁=0. Marks, such as mark 222, indicatewhen the other five sources S2, S3, S4, S5, and S6 are activated withcorresponding dithered time delays τ₂, τ₃, τ₄, τ_(s), and τ₆ after thestart of the shot interval 214. In shot interval 215, not all of thesources are selected for activation and the order in which the sourcesare activated and associated time delays are different from those inshot interval 214. The sources selected for each activation in a shotinterval is called a “source array.” For example, the source array forshot interval 214 is composed of all the sources S1, S2, S3, S4, S5, andS6; and the source array for the shot interval 215 is composed of thesources S1, S2, S3, S4, and S5.

Receivers 118 may include a particle motion sensor that detects particlemotion, velocities, or accelerations over time, a pressure sensor thatdetects variations in water pressure over time, or a combination ofparticle motion and pressure sensors (i.e., a dual sensor). FIG. 3 showsa side-elevation view of the data acquisition system 100 with amagnified view 302 of an example receiver 118. In this example, themagnified view 302 reveals that the receiver 118 is a dual sensorcomposed of a pressure sensor 304 and a particle motion sensor 306. Eachpressure sensor may measure changes in hydrostatic pressure over timeand produces pressure data denoted by p({right arrow over (x)}, t),where {right arrow over (x)} represents the Cartesian coordinates (x, y,z) of the receiver, and t represents time. The particle motion sensorsmay be responsive to water motion. In general, particle motion sensorsdetect particle motion in a direction normal to the orientation of theparticle motion sensor and may be responsive to such directionaldisplacement of the particles, velocity of the particles, oracceleration of the particles. The particle motion sensor data producedby the particle motion sensors may be converted to particle motionvelocity data. For example, when particle motion sensors that areresponsive to position are used, the particle motion sensor data denotedby g_({right arrow over (n)})({right arrow over (x)}, t) may bedifferentiated to convert the data to particle motion velocity datadenoted by v_({right arrow over (n)})({right arrow over (x)}, t), whereunit normal vector {right arrow over (n)} points in the directionparticle motion is measured. Likewise, when particle motion sensors thatare responsive to acceleration (i.e., accelerometers) are used, theparticle acceleration data denoted by a_({right arrow over (n)})({rightarrow over (x)}, t) may be integrated to convert the data to particlemotion velocity data v_({right arrow over (n)})({right arrow over (x)},t). The particle motion sensors are typically oriented so that theparticle motion is measured in the vertical direction (i.e., {rightarrow over (n)}=(0,0,z)) in which case v_(z)({right arrow over (x)}, t)is called the vertical velocity data. Alternatively, each receiver mayinclude two additional particle motion sensors that measure particlemotion in two other directions, {right arrow over (n)}₁ and {right arrowover (n)}₂, that are orthogonal to {right arrow over (n)} (i.e., {rightarrow over (n)}·{right arrow over (n)}₁={right arrow over (n)}·{rightarrow over (n)}₂=0, where “·” is the scalar product) and orthogonal toone another (i.e., {right arrow over (n)}₁·{right arrow over (n)}₂=0).In other words, each receiver may include three particle motion sensorsthat measure particle motion in three orthogonal directions. Forexample, in addition to having a particle motion sensor that measuresparticle motion in the z-direction to give v_(z)({right arrow over (x)},t), each receiver may include a particle motion sensor that measures thewavefield in the in-line direction in order to obtain the inlinevelocity wavefield, t), and a particle motion sensor that measures thewavefield in the cross-line direction in order to obtain the cross-linevelocity wavefield, v_(y)({right arrow over (x)}, t). In certainimplementations, at least some of the receivers may by composed of onlypressure sensors, and in other implementations, at least some of thereceivers may be composed of only particle motion sensors.

The streamers 106-111 and the survey vessels 101-104 may include sensingelectronics and data-processing facilities that allow seismic datagenerated by each receiver to be correlated with the times when thesources are activated, absolute positions on the free surface 112, andabsolute three-dimensional positions with respect to an arbitrarythree-dimensional coordinate system. The pressure data and particlemotion data may be stored at the receiver and/or may be sent along thestreamers and data transmission cables to the survey vessel 101, wherethe data may be stored electronically or magnetically on data-storagedevices located onboard the survey vessel 101. The pressure data andparticle motion data represent pressure and particle motion wavefieldsand, therefore, may also be called the pressure wavefield and particlemotion wavefield, respectively.

In FIG. 3, directional arrow 308 represents the direction of an up-goingwavefield at the location of receiver 310 and dashed-line directionalarrow 312 represents a down-going wavefield produced by an up-goingwavefield reflection from the free surface 112 before reaching thereceiver 310. In other words, the pressure wavefield p({right arrow over(x)}, t) is composed of an up-going pressure wavefield component and adown-going pressure wavefield component, and the particle motionwavefield g_({right arrow over (n)})({right arrow over (x)}, t) iscomposed of an up-going wavefield component and a down-going wavefieldcomponent. The down-going wavefield contaminates pressure and particlemotion data and creates notches in the seismic data spectral domain.Filtering may be done to remove the down-going wavefields from thepressure and particle motion data, leaving the up-going wavefields whichare typically used to analyze the subterranean formation.

The seismic data measured by each pressure sensor or particle motionsensor may be a time series that consist of a number of consecutivelymeasured values called amplitudes separated in time by a sample rate.The time series measured by a pressure or particle motion sensor iscalled a “trace,” which may consist of thousands of samples collected ata typical sample rate of about 1 to 5 ms. A trace is generally arecording of a subterranean formation response to seismic energy thatpasses from an activated source, into the subterranean formation where aportion of the seismic energy is reflected, and ultimately detected by asensor as described above. A trace records variations in atime-dependent amplitude that represents seismic energy in the portionof the wavefield measured by the sensor. In other words, each trace is aset of time-dependent pressure or particle motion sensor amplitudesdenoted by

tr(k)={A _(k)(t _(b))}_(b=1) ^(B)  (1)

where k is a trace index;

A_(k)(t_(b)) is the amplitude of trace k at time sample t_(b); and

B is the number of time samples.

As explained above, the wavefield typically arrives first at thereceivers located closest to the sources. The distance from the sourcesto a receiver is called the “source-receiver offset,” or simply“offset,” which creates a delay in the arrival time of a secondarywavefield from a substantially horizontal interface within thesubterranean formation. A larger offset generally results in a longerdelay in the arrival time. The traces are collected to form a “gather”that can be further processed using various seismic data processingtechniques in order to obtain information about the structure of thesubterranean formation.

FIG. 4 shows example rays that represent paths of an acoustic signal 400output from the source S1. Dashed-line rays, such as rays 402, representseismic energy reflected from the formation surface 122 to the receivers118 located along the streamer 108, and solid-line rays, such as rays404, represent seismic energy reflected from the interface 124 to thereceivers located along the streamer 108. Note that for simplicity ofillustration only a handful of ray paths are represented. Each pressuresensor may measure the hydrostatic pressure and each particle motionsensor may measure the particle motion of the seismic energy reflectedfrom the subterranean formation 120 or interfaces therein. Thehydrostatic pressure data and/or particle motion data generated at eachreceiver may be time sampled and recorded as separate traces. In theexample of FIG. 4, the collection of traces generated by the receiversalong the streamer 108 for a single activation of the source S1 may becollected to form a “common-shot gather” or simply a “shot gather.” Thetraces generated by the receivers located along each of the other fivestreamers for the same activation may be collected to form separate shotgathers, each shot gather associated with one of the streamers.

FIG. 5 shows a plot of a shot gather composed of example traces 501-505of the wavefield measured by the five receivers located along thestreamer 108 shown in FIG. 4. Vertical axis 506 represents time andhorizontal axis 508 represents trace numbers. The traces are arranged sothat trace “1” represents the seismic data generated by the receiverlocated closest to the source S1 and trace “5” represents the seismicdata generated by the receiver located farthest (along the length of thestreamer) from the source S1. The traces 501-505 may represent variationin the amplitude of either the pressure data or the particle motion datameasured by corresponding sensors of the five receivers. The exampletraces include wavelets or pulses 509-513 and 514-518 represented bypeaks colored black that represent the up-going wavefield measured bythe pressure sensors or particle motion sensors. The time intervalsalong the traces 501-505 from the trace number axis 508 (i.e., timezero) to the wavelets 509-513 represents two-way travel time of theseismic energy output from the source S1 to the formation surface 122and to the receivers located along the streamer 108, and the wavelets514-518 represent longer two-way travel time of the seismic energyoutput from the source S1 to the interface 124 and to the same receiverslocated along the streamer 108. The amplitudes of the peaks or troughsof the wavelets 509-513 and 514-518 indicate the magnitude of thereflected seismic energy measured by the receivers.

The arrival times versus source-receiver offset is longer withincreasing source-receiver offset. As a result, the wavelets generatedby a surface or an interface are collectively called a “reflected wave”or simply “reflection” that tracks a parabolic-shaped curve. Forexample, dashed curve 520 represents the distribution of the wavelets510-513 reflected from the formation surface 122, which are called a“surface reflected wave” or “water bottom reflection” and solidparabolic curve 522 represents the distribution of the wavelets 511-515from the interface 124, which are called an “interface reflected wave”or “interface reflection.”

FIG. 6 shows an example of a gather composed of 38 traces. Each trace,such as trace 602, varies in amplitude over time and represents seismicenergy reflected from the surface and five different interfaces within asubterranean formation as measured by a pressure sensor or a particlemotion sensor. In the expanded view, wavelets that correspond toreflection from the same surface or interface of the subterraneanformation appear chained together. For example, wavelets 604 with theshortest transit time represent a water bottom reflection, and wavelets606 represent an interface reflected wave emanating from an interfacejust below the surface. Reflected waves 608-611 represent reflectionsfrom interfaces located deeper within the subterranean formation.

The gathers shown in FIGS. 5 and 6 are described for seismic data sortedinto a common-shot domain (or simply the shot domain). A domain is acollection of gathers that share a common attribute with respect to theseismic data recording locations. The seismic data may be sorted intoany suitable domain for examining the features of a subterraneanformation including, for example, a common-receiver domain,common-receiver-station domain, or common-midpoint (“CMP”) domain.

FIG. 7 shows a plot of different ways seismic data collected in a surveymay be sorted into domains. Vertical axis 702 represents the in-linereceiver coordinates and horizontal axis 704 represents the in-linesource coordinates. X's, such as X 706, represent where a recording(i.e., pressure or particle motion) has taken place. In this plot, acolumn of recordings identified by dashed line 708 represents a shotgather, and a row of recordings identified by dashed line 710 representsa common-receiver-station gather. Recordings collected along a diagonalrepresented by dashed line 712 is a common-receiver gather, andrecordings collected along a diagonal represented by dashed line 714 isa CMP gather. The gathers form different domains. For example, the shotgathers form a shot domain, the common-receiver gathers form acommon-receiver domain, the common-receiver-station gathers form acommon-receiver-station domain, and the CMP gathers form a CMP domain.Certain domains are orthogonal. For example, as shown in FIG. 7, thegathers in the shot domain are orthogonal to the gathers in thecommon-receiver domain.

Seismic data processing may include correction of the effects ofdifferent source-receiver offsets using a process called “normalmoveout” (“NMO”) or a process called “linear moveout” (“LMO”) andcorrection of angled interfaces in a process called “dip moveout”(“DMO”). A dip is an angle of inclination with respect to a horizontalplane. For example, in FIG. 4 angle θ 406 represents a dip angle belowhorizontal. NMO and DMO corrections have the effect of flattening outthe parabolic-shaped reflections and may be applied as follows. Prior toNMO and DMO corrections NMO velocity analysis may be performed on agather, such as a CMP gather, to determine a regional velocity model. Afirst NMO correction may be applied to the data using the regionalvelocity model followed by application of a DMO correction. An inverseNMO correction may be applied using the regional velocity functionfollowed by performing a more detailed NMO velocity analysis todetermine a detailed velocity model. Finally, NMO correction may beapplied again using the detailed velocity model to give a gather withflatten reflections. After NMO and DMO corrections have been performedto flatten the reflections, traces from different shot records with acommon reflection point may be stacked to form a single trace duringseismic data processing. For example, NMO and DMO corrections may beapplied to CMP gathers followed by stacking in order to improve thesignal-to-noise ratio, reduce noise, improve seismic data quality, andreduce the amount of data. LMO may be used with an averagesource-to-receiver velocity model.

Dithered activation of two or more seismic sources as described abovewith reference to FIG. 2 allows for the attenuation of the seismicenergy from interfering sources after sorting the acquired seismic datato an appropriate domain, such as the domains described above withreference to FIG. 7. After sorting to an appropriate domain, the seismicdata associated with one of the sources (i.e., the primary source) maybe aligned with time zero and appears coherent, while the seismic dataassociated with the one or more other sources (i.e., secondary sources)appears incoherent. As a result, methods that discriminate betweencoherent and incoherent energy allow for the sources to be separated.

As explained above, blended seismic data is composed of interferingseismic energy emitted from a subterranean formation for each of the twoor more source activations. Methods and systems separate, or “deblend,”the blended seismic data into separate deblended seismic data gathers(i.e., deblended gathers) that are each associated with activation ofone of the two or more sources as if each deblended gather is generatedwithout substantial interference from seismic energy generated by theother sources. For example, a frequency-domain mathematical relationshipbetween the blended seismic data and the deblended gathers may berepresented as follows:

P′(z _(r) ,z _(s))=P(z _(r) ,z _(s))Γ  (2)

where P′(z_(r), z_(s)) is a frequency-domain blended seismic datamatrix;

-   -   P(z_(r), z_(s)) is a frequency-domain deblended seismic data        matrix;    -   r is a frequency-domain blending matrix; and    -   z_(r) and z_(s) are the receiver and source depth levels,        respectively.

FIG. 8 shows a more detailed mathematical representation of the matricesof Equation (2). Matrix 801 represents the matrix P′(z_(r), z_(s)) ofblended seismic data, matrix 802 represents matrix P(z_(r),z_(s)) ofdeblended seismic data, and matrix 803 represents blending matrix Γ.Blended seismic data matrix P′(z_(r),z_(s)) 801 is an R×L matrix,deblended seismic data matrix P(z_(r), z_(s)) 802 is an R×S matrix, andblending matrix Γ803 is an S×L matrix, where R represents the totalnumber of receivers in the data acquisition surface, S is the totalnumber of sources, and L is the number of sources arrays activated in asurvey. Consider elements P_(rs) of the deblended seismic data matrixP(z_(r), z_(s)) 802. Each element P_(rs) represents a deblended trace,where subscript r is an integer receiver index 1≦r≦R, and subscript s isan integer source index 1≦s≦S. Each row in the deblended seismic datamatrix P(z_(r),z_(s)) 802 corresponds to a different receiver and eachcolumn corresponds to a different source. For example, matrix elementsP₁₁, P₁₂, and P₁₃ in row 804 represent traces generated by receiver 1 asa result of separate, or deblended, seismic energy generated bycorresponding sources S1, S2, and S3. Matrix elements P₁₁, P₂₁, and P₃₁in column 805 represent traces generated by receivers 1, 2, and 3 as aresult of separate, or deblended, seismic energy generated by the sourceS1. Consider matrix elements P_(si) of the blending matrix Γ 803, wheresubscript l is an integer source array index 1≦l≦L. Because eachrecording period has a different source array, the subscript l may beused as a recording period index. For surveys with dithered activationtimes and sources with substantially equal seismic energy output, thematrix elements of the blending matrix Γ 803 are time-delay phase givenby δ_(sl)=exp(−jωτ_(sl)), where τ_(sl) is the time delay for source s insource array l. Each row in the blending matrix Γ 803 corresponds to adifferent source and each column corresponds to a different sourcearray. Finally, consider matrix elements P′_(rl) in blended seismic datamatrix P′(z_(r), z_(s)) 801. Each element P′_(rl) represents a tracegenerated by a receiver r as a result of seismic energy produced by asource array l. Each column in the blended seismic data matrixP′(z_(r),z_(s)) 801 is composed of R blended traces generated by all Rreceivers of the data acquisition surface in a recording period as aresult of activating a particular source array. Traces in the blendedseismic data matrix P′(z_(r), z_(s)) 801 may be sorted into differentdomains by selecting a particular set of traces that correspond to thedomain. For example, common-shot gathers are created by selecting tracesin the same column; common-receiver gathers are created be selectedtraces along a diagonal.

FIG. 9 shows a specific example of a matrix equation that may beobtained using the data acquisition system 100 in a single recordingperiod. As shown in FIG. 1B, the example data acquisition system 100contains R=30 receivers and S=6 sources. Returning to FIG. 9, blendingmatrix 901 is obtained from the time delays determined after thebeginning of a shot interval. For example, the source array for the shotinterval 215 in FIG. 2 is composed of sources S1, S2, S3, S4, and S5with associated time delays τ₁₁, τ₂₁, τ₃₁, τ₄₁=0 (i.e., the shotinterval starts with activating source S4), and τ₅₁. The correspondingfrequency-domain blending matrix 901 is given by

$\begin{matrix}{\Gamma = {\begin{bmatrix}\Gamma_{11} \\\Gamma_{21} \\\Gamma_{31} \\\Gamma_{41} \\\Gamma_{51} \\\Gamma_{61}\end{bmatrix} = \begin{bmatrix}{\exp \left( {{- j}\; \omega \; \tau_{11}} \right)} \\{\exp \left( {{- j}\; \omega \; \tau_{21}} \right)} \\{\exp \left( {{- j}\; \omega \; \tau_{31}} \right)} \\1 \\{\exp \left( {{- j}\; \omega \; \tau_{51}} \right)} \\0\end{bmatrix}}} & (3)\end{matrix}$

Note that because the source S6 is not activated in the source array,the corresponding matrix element Γ₆₁ is equal to zero. Blended seismicdata matrix 902 is composed of 30 traces generated by all 30 receiversas a result of activating an array of sources with time delaysrepresented in the blending matrix 901. Elements P_(rs) of the deblendedseismic data matrix 903 represent recorded seismic energy generated by asource s and measured at a receiver r without interference of seismicenergy generated by the other five sources.

Returning to FIG. 8, the traces in the blended seismic data matrixP′(z_(r), z_(s)) 801 are obtained from recording seismic data at thereceivers and elements of the blending matrix Γ 803 are determined fromthe known times when the sources are activated in a shot interval. Inthe frequency domain, Equation (2) indicates that each blended traceelement in the blended seismic data matrix P′(z_(r), z_(s)) 801 is alinear combination of a row of deblended traces in deblended seismicdata matrix P(z_(r), z_(s)) 802 and a column of phase terms in theblending matrix Γ 803:

$\begin{matrix}{P_{rl}^{\prime} = {\sum\limits_{s = 1}^{S}\; {P_{rs}\Gamma_{sl}}}} & (4)\end{matrix}$

However, the traces P_(rs) in the deblended seismic data matrixP(z_(r),z_(s)) 802 are unknown and may not be measured in isolation.Ideally, traces in deblended seismic data matrix P(z_(r),z_(s)) 802 maybe determined from matrix inversion:

P(z _(r) ,z _(s))=P′(z _(r) ,z _(s))Γ⁻¹  (5)

where Γ⁻¹ is an inverse blending matrix of the blending matrix Γ.

However, in practice, it may be the case that the matrix equation inEquation (2) cannot be solved with matrix inversion according toEquation (5) because the matrix equation in Equation (2) isunderdetermined (i.e., L<R). In cases where L<R, the blending matrix Γis not invertible and a unique deblended seismic data matrix P(z_(r),z_(s)) 802 does not exist.

Because the deblended seismic data matrix P may not be directlydetermined from matrix inversion, methods of the current applicationiteratively determine a deblended data matrix M, which is anapproximation or model of the deblended seismic data matrix P. Thedeblended data matrix M may be determined from a residual equation thatcharacterizes the closeness of a blended seismic data matrix P′ to thedeblended data matrix M with time delays represented by blending matrixΓ:

RES=P′−MΓ  (6)

An example of a deblended data matrix M is an R×S matrix illustrated inFIG. 10 and is an approximation of the deblended seismic data matrix P.Each element in the deblended data matrix M is a deblended trace denotedby M_(rs). The columns of the deblended data matrix M are denoted bym_(s) and each column is a deblended gather associate with activation ofone source s. Each deblended gather m_(s) is an approximation or modelof seismic data measured for one source of a source array activated in ashot interval without interfering seismic energy created by the othersources in the source array. For example, rectangles 1001-1004 representS deblended gathers m₁, m₂, m₃, and m_(s). Each deblended gather m_(s)is a column of R traces in the deblended data matrix M. For example, asshown in FIG. 10, deblended gather m_(s) 1004 is composed of the columnof traces M_(rs), where 1≦r≦R.

The magnitude of the residual RES may be iteratively minimized whilealso iteratively constructing the deblended data matrix M by redefiningthe residual in Equation (6) as an iterative equation:

RES _(i) =P′−M _(i)Γ  (7)

where subscript i is a non-negative integer iteration index; and

M_(i)=[m_(1,i) m_(2,i) m_(3,i) . . . m_(s,i)] with each m_(s,i) is adeblended gather associated with activation of a source s.

In the time domain, the iterative equation represented by Equation (7)becomes:

res _(i) ={circumflex over (P)}′−{circumflex over (M)} _(i)*{circumflexover (Γ)}  (8)

where {circumflex over (P)}′ is a time-domain blended seismic datamatrix;

-   -   {circumflex over (M)}_(i) is a time-domain deblended data matrix        characterized by {circumflex over (M)}_(i)=[{circumflex over        (m)}_(1,i) {circumflex over (m)}_(2,i) {circumflex over        (m)}_(3,i) . . . {circumflex over (m)}_(S,i)] with each        {circumflex over (m)}_(s,i) a time-domain deblended gather        associated with activation of a source s;    -   {circumflex over (Γ)} is a time-domain blending matrix;    -   res_(i) is a time-domain residual; and    -   * represents convolution.

FIGS. 11 and 12 present an example method for separating blended seismicdata {circumflex over (P)}′ obtained using SSA in the time domain. FIG.11 shows an example flow-control diagram of an iterative method forseparating seismic data obtained using SSA. In block 1101, blendedseismic data {circumflex over (P)}′ may be received, wherein blendedseismic data {circumflex over (P)}′ represents energy measured byreceivers of a data acquisition surface in L recording periods. Block1101 also includes receiving the time delays τ_(sl) associated withactivation of a source array in each shot interval of the L recordingperiods.

In block 1102, parameters are initialized. The iteration index i isinitialized to zero (i.e., i=0) and the residual for i=0 is initializedto the blended seismic data matrix:

res ₀ ={circumflex over (P)}′  (9)

and the deblended data matrix is initialized to

M ₀ =[{circumflex over (m)} _(1,0) {circumflex over (m)} _(2,0){circumflex over (m)} _(3,0) . . . {circumflex over (m)} _(S,0)]=[0 0 0. . . 0]  (10)

which represents the amplitudes of each deblended trace in the initialdeblended data matrix {circumflex over (M)}₀ are initialized to zero(i.e., {circumflex over (M)}_(rs)=0, for 1≦r≦R and 1≦s≦S). Parameterinitialization also includes initializing a semblance-value threshold,T, to a value in the interval 0<T≦1. The semblance-value threshold T maybe input to a routine called in block 1104. For example, the thresholdmay be initialized to 0.6. In certain implementations, the threshold Tmay be fixed for each iteration. In other implementations, the thresholdmay be changed after each iteration. For the particular implementationillustrated in FIG. 11, the operations executed in blocks 1103-1106include changing the threshold in block 1106 as explained below.

In block 1103, the residual res_(i) may be sorted into a particulardomain, such as any one of the domains described above with reference toFIG. 7. For example, for each iteration i, the residual res_(i) may besorted into the common-receiver domain. In other implementations, theresidual res_(i) may be sorted into different domains for certainiterations. For example, for the first iteration with i=0 the residualres₀ may be sorted into the common-receiver domain, and in a subsequentiteration with i=1, the residual res₁ may be sorted into the CMP domain.

In block 1104, a subroutine “update residual and deblended seismic datagathers” may be called to receive as input the residual and the set ofdeblended gathers and output an updated residual and an updated set ofdeblended gathers as described below with reference to FIG. 12. Indecision block 1105, when the iteration index i is greater than a userselected number of iterations, N, control flows to block 1107 in whichthe most recent set of deblended gathers {{circumflex over(m)}_(s,i)}_(s=1) ^(S) is output. Each gather m_(s,i) output in block1107 may be substantially free of seismic energy contamination fromother sources. Otherwise, control flows to block 1106. In block 1106,the semblance-value threshold T is decreased. For example, the thresholdT may be initialized to a value 0.6 for iteration i=0 decreased to avalue 0.4 for iteration i=1 and further decreased to a value 0.2 foriteration i=2.

Implementations are not limited to using a fixed number of iterations Nto determine when the iteration process in blocks 1103-1106 ends indecision block 1105. In other implementations, the decision in decisionblock 1105 to end the iterative process in blocks 1103-1106 may be basedon when the magnitude of the difference between a current residualres_(i) and a previous residual res_(i-1) is less than a selectedresidual threshold, τ_(res). For example, the decision block 1105 may beimplemented as follows. When

τ_(r) ≧∥res _(i) −res _(i-1)∥  (11)

control flows to block 1107, otherwise, control flows to block 1106.

FIG. 12 shows a flow-control diagram for the routine “update residualand deblended seismic data” called in block 1104 of the flow-controldiagram of FIG. 11. In block 1201, the iteration index i is incremented.In block 1202, pseudo-deblending of the residual res_(i-1) may beperformed for each time delay τ_(si) by subtracting the time delay fromeach trace of the residual res_(i-1) to give a correspondingpseudo-deblended residual res_(i-1) ^((s)). Pseudo-deblending theresidual may be performed for each activated source to give a set ofpseudo-deblended residuals:

res _(i-1) ⁽¹⁾ res _(i-1) ⁽²⁾ res _(i-1) ⁽³⁾ , . . . , res _(i-1)^((s))  (12)

Pseudo-deblending creates pseudo-deblended residuals res_(i-1) ^((s))that are each aligned in time with a time delay τ_(sl). In other words,each pseudo-deblended residual res_(i-1) ^((s)) is aligned in time witha gather m_(s,i-1). In block 1203, velocity analysis and NMO or LMO andDMO corrections may be applied as described above. In particular, whenthe seismic data is sorted into the CMP domain in block 1103 of FIG. 11,velocity analysis, NMO, and DMO may be applied to flatten reflectors.

The diagram in FIG. 12 includes an outer for-loop beginning in block1204 repeats the operations represented by blocks 1205-1212 for each ofthe S sources activated in the shot interval. An inner for-loopbeginning in block 1205 repeats the operations represented by blocks1206-1209 for each time sample t_(b) described above with reference toEquation (2).

In block 1206, semblance values S_(v)(t_(b)) may be calculated for eachtime sample level t_(b). For example, let tr denote the number of tracesin the pseudo-deblended residual res_(i-1) ^((s)), and let a_(k)(t_(b))be the amplitude of trace k at time sample level t_(b) inpseudo-deblended residual res_(i-1) ^((s)). In one implementation, thesemblance value S_(v)(t_(b)) may be calculated using overlapping windowsthat are Q traces wide for each time sample level t_(b) according to thefollowing pseudo-code:

$\begin{matrix}{{{{set}\mspace{14mu} {window}\mspace{14mu} {width}\mspace{14mu} Q\mspace{14mu} {to}\mspace{14mu} {an}\mspace{14mu} {odd}\mspace{14mu} {number}\mspace{14mu} {of}\mspace{14mu} {traces}\mspace{14mu} \left( {Q < {tr}} \right)};}{{{for}\mspace{14mu} {each}\mspace{14mu} {time}\mspace{14mu} {sample}\mspace{14mu} {level}\mspace{14mu} t_{b}};}\mspace{20mu} {{{{for}\mspace{14mu} {each}\mspace{14mu} q} = 1},\ldots \mspace{14mu},{{tr};}}\mspace{40mu} {{S_{v}\left( t_{b} \right)} = \frac{\left( {\sum\limits_{k = {q - {{({Q - 1})}/2}}}^{q + {{({Q - 1})}/2}}\; {a_{k}\left( t_{b} \right)}} \right)^{2}}{{Q\left( {\sum\limits_{k = {q - {{({Q - 1})}/2}}}^{q + {{({Q - 1})}/2}}\; {a_{k}\left( t_{b} \right)}} \right)}^{2}}}} & (13)\end{matrix}$

In block 1207, when semblance value S_(v)(t_(b)) is greater than orequal to the semblance-value threshold T, control flows to block 1208;otherwise, control flows to decision block 1209. The coherency value isa measure of coherence in the amplitudes a_(k)(t_(b)) of thepseudo-deblended residual res_(i-1) ^((s)) at time sample t_(b). Inblock 1208, a coherency filter may be applied to the set of amplitudes{a_(k)(t_(b))}, where k ranges from

$q - {\frac{\left( {Q - 1} \right)}{2}\mspace{14mu} {to}\mspace{14mu} q} + {\frac{\left( {Q - 1} \right)}{2}.}$

The coherency filter may be any one of a median filter, frequency-spacedeconvolution, tau-p filter, or a frequency-wavenumber filter used toseparate coherent energy from incoherent energy. In decision block 1209,when all of the time samples of pseudo-deblended residual res_(i-1)^((s)) have been considered, control flows to block 1210; otherwise theoperations represented blocks 1206-1208 are repeated for another timesample level.

In block 1210, the coherent energy obtained in block 1208 for all timesamples of the pseudo-deblended residual res_(i-1) ^((s)) may becollected to form a corresponding coherent-energy gather, c^((s)). Inblock 1211, the coherent energy of coherent-energy gather c^((s)) may beused to update a corresponding previous deblended gather to give anupdated deblended gather:

{circumflex over (m)} _(s,i) ={circumflex over (m)} _(s,i-1) +c^((s))  (14)

In decision block 1212, when the operations in blocks 1205-1212 havebeen repeated for all S, control flows to block 1214. Otherwise, controlflows to block 1213 in which s is incremented and the operations inblocks 1205-1212 are repeated.

In block 1214, the residual may be updated by first blending thedeblended seismic data to obtain blended coherent energy:

$\begin{matrix}{{{blend}\left( {c^{(1)},c^{(2)},\ldots \mspace{14mu},c^{(S)}} \right)} = {\sum\limits_{s = 1}^{S}\; {{time\_ shift}_{\tau_{sl}}\left( c^{(s)} \right)}}} & (15)\end{matrix}$

where time_shift(c^((s))) is calculated by adding τ_(sl) to the timecomponent of each trace in the coherent-energy gather c^((s)).

The blended coherent energy blend (c⁽¹⁾, c⁽²⁾, . . . , c^((s))) may besubtracted from the residual to give an updated residual:

res _(i) =res _(i-1) ^((s))−blend(c ⁽¹⁾ ,c ⁽²⁾ , . . . , c ^((s)))  (16)

After the residual res_(i) has been updated, the method returns toexecute the operation represented by decision block 1105 in theflow-control diagram in FIG. 11.

FIG. 13 shows an example of a generalized computer system that executesefficient methods for seismic source separation and therefore representsa geophysical-analysis data-processing system. The internal componentsof many small, mid-sized, and large computer systems as well asspecialized processor-based storage systems can be described withrespect to this generalized architecture, although each particularsystem may feature many additional components, subsystems, and similar,parallel systems with architectures similar to this generalizedarchitecture. The computer system contains one or multiple centralprocessing units (“CPUs”) 1302-1305, one or more electronic memories1308 interconnected with the CPUs by a CPU/memory-subsystem bus 1310 ormultiple busses, a first bridge 1312 that interconnects theCPU/memory-subsystem bus 1310 with additional busses 1314 and 1316, orother types of high-speed interconnection media, including multiple,high-speed serial interconnects. The busses or serial interconnections,in turn, connect the CPUs and memory with specialized processors, suchas a graphics processor 1318, and with one or more additional bridges1320, which are interconnected with high-speed serial links or withmultiple controllers 1322-1327, such as controller 1327, that provideaccess to various different types of computer-readable media, such ascomputer-readable medium 1328, electronic displays, input devices, andother such components, subcomponents, and computational resources. Theelectronic displays, including visual display screen, audio speakers,and other output interfaces, and the input devices, including mice,keyboards, touch screens, and other such input interfaces, togetherconstitute input and output interfaces that allow the computer system tointeract with human users. Computer-readable medium 1328 is adata-storage device, including electronic memory, optical or magneticdisk drive, USB drive, flash memory and other such data-storage device.The computer-readable medium 1328 can be used to store machine-readableinstructions that encode the computational methods described above andcan be used to store encoded data, during store operations, and fromwhich encoded data can be retrieved, during read operations, by computersystems, data-storage systems, and peripheral devices.

FIGS. 14A-14C show results of seismic source separation applied toblended seismic data generated by two sources sorted into common-shotgathers. FIG. 14A shows a common-shot gather of blended seismic dataobtained from SSA using a primary source and a secondary source. Theenergy generated by the secondary source is identified by dashed oval1401. FIG. 14B shows the common-shot gather after applying typicalseismic source separation, which does not include a semblance valueconstraint, to give deblended seismic data associated with the primarysource. However, dashed oval 1402 encircles a significant amount ofinterfering energy originating from the secondary source that has leakedpast the filter to contaminate the energy associated with the primarysource. FIG. 14C shows the common-receiver gather after applying theseismic source separation method with semblance analysis described aboveto give deblended seismic data associated with the primary source.Dashed oval 1403 encircles a small portion of interfering energy thatoriginates from the secondary source and has leaked passed the filter.The results displayed in FIG. 14C reveal that the method described aboveallows significantly less energy originating from the secondary sourceto contaminate the energy associated with the primary source.

The methods described above may be implemented in real time on board asurvey vessel while a survey is being conducted. For example, deblendedgathers may be generated after each recording period. The deblendedgathers described above may produce a geophysical data productindicative of certain properties of a subterranean formation. Thegeophysical data product may include processed seismic geophysical dataand may be stored on a computer-readable medium as described above. Thegeophysical data product may be produced offshore (i.e. by equipment onsurvey vessel) or onshore (i.e. at a computing facility on land) eitherwithin the United States or in another country. When the geophysicaldata product is produced offshore or in another country, it may beimported onshore to a data-storage facility in the United States. Onceonshore in the United States, geophysical analysis may be performed onthe data product.

Although the above disclosure has been described in terms of particularembodiments, it is not intended that the disclosure be limited to theseembodiments. Modifications within the spirit of the disclosure will beapparent to those skilled in the art. For example, any of a variety ofdifferent implementations of the method described above for seismicsource separation may be carried out by varying any of many differentdesign and development parameters, including programming language,underlying operating system, modular organization, control structures,data structures, and other such design and development parameters. Themethods may also be carried out in the frequency domain by firsttransforming the blended seismic data {circumflex over (P)}′ in the timedomain to blended seismic data {circumflex over (P)}′ in the frequencydomain, using multiplication by the time-delay phaseΓ_(si)=exp(−jωτ_(sl)) to affect time shifts, and transforming theresulting frequency-domain deblended gathers {m_(s,i)}_(s=1) ^(S) backto time-domain deblended gathers {{circumflex over (m)}_(s,i)}_(s=1)^(S). Examples of suitable transforms include, but are not limited to, afast Fourier transform and a Laplace transform.

It is appreciated that the previous description of the disclosedembodiments is provided to enable any person skilled in the art to makeor use the present disclosure. Various modifications to theseembodiments will be readily apparent to those skilled in the art, andthe generic principles defined herein may be applied to otherembodiments without departing from the spirit or scope of thedisclosure. Thus, the present disclosure is not intended to be limitedto the embodiments shown herein but is to be accorded the widest scopeconsistent with the principles and novel features disclosed herein.

1. A method for deblending seismic data generated by multiple seismicsources, the method comprising: receiving blended seismic data recordedby receivers, the blended seismic data generated by simultaneousactivation of two or more seismic sources with dithered time delays, thetwo or more sources located at different positions; calculating aresidual based on the blended seismic data and the dithered time delays;calculating two or more pseudo-deblended residuals, eachpseudo-deblended residual is the residual aligned in time with anactivation time of one of the two or more sources; and iterativelycalculating deblended gathers from the pseudo-deblended residuals, eachdeblended gather represents seismic energy produced by one of the two ofmore sources.
 2. The method of claim 1 executed on a programmablecomputer programmed to perform the method.
 3. The method of claim 1wherein calculating the residual further comprises calculating adifference between the blended seismic data and deblended gathersshifted by corresponding time-delay.
 4. The method of claim 1 whereincalculating two or more pseudo-deblended residuals further comprisecalculating each pseudo-deblended residual from the residual shifted byone of the time delays.
 5. The method of claim 1 wherein iterativelycalculating the deblended gathers from the pseudo-deblended residualsfurther comprises: initializing amplitudes of deblended gathers to zero;for each pseudo-deblended residual, calculating a semblance value foramplitudes of the pseudo-deblended residual at each time-sample level;applying a coherence filter to amplitudes with a semblance value greaterthan a semblance-value threshold at each time-sample level; forming acoherent-energy gather from the coherent energy filtered from thepseudo-deblended residual by the coherency filter; and updating adeblended gather associated with the pseudo-deblended residual by addingthe coherent-energy gather to the deblended gather.
 6. The method ofclaim 5 further comprises decreasing the semblance-value threshold foreach calculation of the deblended gathers.
 7. The method of claim 5further comprises: blending the coherent-energy gathers obtained foreach of the pseudo-deblended residuals to generate a blendedcoherent-energy gather; and subtracting the blended coherent-energygather from the residual to generate an updated residual.
 8. The methodof claim 1 further comprises sorting the blended seismic data into adomain.
 9. The method of claim 1 further comprises storing the deblendedgathers in one or more data-storage devices.
 10. The method of claim 1further comprising producing a geophysical data product from thedeblended gathers.
 11. The method of claim 10, further comprisingrecording the geophysical data product on a physical, non-volatilecomputer-readable medium suitable for importing onshore.
 12. The methodof claim 11, further comprising performing geophysical analysis onshoreon the geophysical data product.
 13. A computer system for deblendingseismic data generated by multiple seismic sources, the systemcomprising: one or more processors; one or more data-storage devices;and a routine stored in one or more of data-storage devices and executedby the one or more processors, the routine directed to receiving blendedseismic data recorded by receivers, the blended seismic data generatedby simultaneous activation of two or more seismic sources with ditheredtime delays, the two or more sources located at different positions;calculating a residual based on the blended seismic data and thedithered time delays; calculating two or more pseudo-deblendedresiduals, each pseudo-deblended residual is the residual aligned intime with an activation time of one of the two or more sources; anditeratively calculating deblended gathers from the pseudo-deblendedresiduals, each deblended gather represents seismic energy produced byone of the two of more sources.
 14. The computer system of claim 13wherein calculating the residual further comprises calculating adifference between the blended seismic data and deblended gathersshifted by corresponding time-delay.
 15. The computer system of claim 13wherein calculating two or more pseudo-deblended residuals furthercomprise calculating each pseudo-deblended residual from the residualshifted by one of the time delays.
 16. The computer system of claim 13wherein iteratively calculating the deblended gathers from thepseudo-deblended residuals further comprises: initializing amplitudes ofdeblended gathers to zero; for each pseudo-deblended residual,calculating a semblance value for amplitudes of the pseudo-deblendedresidual at each time-sample level; applying a coherence filter toamplitudes with a semblance value greater than a semblance-valuethreshold at each time-sample level; forming a coherent-energy gatherfrom the coherent energy filtered from the pseudo-deblended residual bythe coherency filter; and updating a deblended gather associated withthe pseudo-deblended residual by adding the coherent-energy gather tothe deblended gather.
 17. The computer system of claim 16 furthercomprises decreasing the semblance-value threshold for each calculationof the deblended gathers.
 18. The computer system of claim 16 furthercomprises: blending the coherent-energy gathers obtained for each of thepseudo-deblended residuals to generate a blended coherent-energy gather;and subtracting the blended coherent-energy gather from the residual togenerate an updated residual.
 19. The computer system of claim 13further comprising producing a geophysical data product from thedeblended gathers.
 20. The computer system of claim 19, furthercomprising recording the geophysical data product on a physical,non-volatile computer-readable medium suitable for importing onshore.21. The computer system of claim 20, further comprising performinggeophysical analysis onshore on the geophysical data product.
 22. Aphysical computer-readable medium having machine-readable instructionsencoded thereon for enabling one or more processors of a computer systemto perform the operations of receiving blended seismic data recorded byreceivers, the blended seismic data generated by simultaneous activationof two or more seismic sources with dithered time delays, the two ormore sources located at different positions; calculating a residualbased on the blended seismic data and the dithered time delays;calculating two or more pseudo-deblended residuals, eachpseudo-deblended residual is the residual aligned in time with anactivation time of one of the two or more sources; and iterativelycalculating deblended gathers from the pseudo-deblended residuals, eachdeblended gather represents seismic energy produced by one of the two ofmore sources.
 23. The computer-readable medium of claim 22 whereincalculating the residual further comprises calculating a differencebetween the blended seismic data and deblended gathers shifted bycorresponding time-delay.
 24. The computer-readable medium of claim 22wherein calculating two or more pseudo-deblended residuals furthercomprise calculating each pseudo-deblended residual from the residualshifted by one of the time delays.
 25. The computer-readable medium ofclaim 22 wherein iteratively calculating the deblended gathers from thepseudo-deblended residuals further comprises: initializing amplitudes ofdeblended gathers to zero; for each pseudo-deblended residual,calculating a semblance value for amplitudes of the pseudo-deblendedresidual at each time-sample level; applying a coherence filter toamplitudes with a semblance value greater than a semblance-valuethreshold at each time-sample level; forming a coherent-energy gatherfrom the coherent energy filtered from the pseudo-deblended residual bythe coherency filter; and updating a deblended gather associated withthe pseudo-deblended residual by adding the coherent-energy gather tothe deblended gather.
 26. The computer-readable medium of claim 25further comprises decreasing the semblance-value threshold for eachcalculation of the deblended gathers.
 27. The computer-readable mediumof claim 25 further comprises: blending the coherent-energy gathersobtained for each of the pseudo-deblended residuals to generate ablended coherent-energy gather; and subtracting the blendedcoherent-energy gather from the residual to generate an updatedresidual.
 28. The computer-readable medium of claim 25 furthercomprising producing a geophysical data product from the deblendedgathers.
 29. The computer-readable medium of claim 28, furthercomprising recording the geophysical data product on a physical,non-volatile computer-readable medium suitable for importing onshore.30. The computer-readable medium of claim 29, further comprisingperforming geophysical analysis onshore on the geophysical data product.